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The dynamics of thermally excited vortices in a dilute two-dimensional Josephson junction array 
where a fraction of the superconducting islands is missing has been investigated using a multiple 
trapping model. An expression for the frequency dependent mobility of vortices has been calculated 
which allows to obtain the frequency dependent complex electrodynamic response of the array for 
different fractions of missing islands. 

PACS numbers: 74.50.+r, 74.60.Ge, 05.90.+m 
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I. INTRODUCTION 

Josephson junction arrays (JJA) consist of supercon- 
ducting islands which are usually arranged on an or- 
dered lattice and coupled by Josephson interaction. Two- 
dimensional (2d) arrays offer a unique opportunity for 
studying a variety of topics in 2d physics, such as phase 
transitions, non-linear dynamics, percolation, frustration 
and disorder, in relatively clean experimental realisation. 
Fabrication of arrays and their basic physical properties 
have been described in various articles. 1,2 The islands be- 
come superconducting at a given transition temperature 
T°. Below this temperature each island I is characterized 
by its superconducting wave function 



ipl = \ipi\e 1 ' 



(1) 



with its amplitude and its phase. For all matters and 
purposes one can assume that \tpi\ has the same value 
in each island, such that the phase is the only relevant 
variable. The islands are linked to each other by the 
Josephson coupling. The potential energy of the array is 
then given by 



H = J2 J »' 






cos(fy - v v 



(2) 



(W) 



The sum can usually be restricted to nearest neighbors 
in the array and the corresponding Josephson coupling J 
is related to the critical current I c by 



2e 



(3) 



A charging energy required when Cooper pairs move 
from one island to another has to be added to the expres- 
sion (2). It is given by a capacitance matrix Cu> coupling 
the time derivatives of the phases (respectively their con- 
jugate momenta) on different sites. 1,2 ' 3 Arrays can then 
be divided into classical and quantum arrays depending 
on the ratio of the Josephson coupling energy to the rel- 
evant charging energy. The experimental work we are re- 
ferring to in this article has been performed on classical 
arrays for which charging energies are unimportant. Dy- 
namic phenomena for such arrays are usually described in 
the framework of the resistively shunted junction (RSJ) 



model. 2,3 ' 4 A resistance matrix Rw describes the normal 
currents flowing between the islands. 

Research on disordered JJAs is a particularly exciting 
field. In an array where some sites are missing the elec- 
trical and magnetic properties of the lattice may change 
drastically depending on the number of missing sites rel- 
ative to the whole array. This effect is related to per- 
colation. About half a century ago this concept was 
introduced by Flory 5 and Stockmayer 6 in order to de- 
scribe disordered structures. Later on the geometrical 
and statistical concepts introduced by Broadbend and 
Hammersley in connection with the study of diffusion of 
fluids in random media, introduction of fractal concepts 8 
and the development of computers and simulation tech- 
niques as well as theoretical and experimental analysis 
have contributed to deepening our understanding of the 
phenomenon of percolation. 

We consider an infinite triangular or square lattice, 
where each site is randomly occupied by a superconduct- 
ing island with some probability p and empty with prob- 
ability 1 — p. The electrical and magnetic properties of 
the array depend on the value of this percolation frac- 
tion p. At a certain value p c of p the superconducting 
properties of the array are destroyed. p c is called crit- 
ical percolation limit. At high concentration of missing 
sites, the occupied sites are either isolated or form small 
clusters of nearest-neighbors, and no path connecting op- 
posite edges of the lattice exists. With increasing p, the 
mean cluster size of the occupied sites increases up to the 
appearance of an infinite cluster which connects opposite 
edges of the whole lattice. The existence of an infinite 
percolating cluster of occupied sites is connected with the 
criticality of the percolation (p c )- When p is further in- 
creased, more and more sites become part of the infinite 
cluster of occupied sites, and finally at p — 1, all sites 
belong to the infinite cluster. The percolation threshold 
depends on different aspects of the lattice structure like 
its form (triangular, square or others), dimensionality, 
and on the type of percolation (bond, site, etc.). For ex- 
ample, in a triangular 2d array, for site percolation, one 
can analytically 9 show that p c = 0.5, but this value can 
be different for other dimensions of the array and other 
types of disorder. 

Different types of disorder can be considered: 



a. For bond disorder the parameters characterizing 
the individual junctions, namely the Josphson couplings 
Jw, the junction resistances Rw and the mutual capac- 
ities Cw , are random. This type of disorder is realized 
when the positions of the individual superconducting is- 
lands deviate in a random way from their ideal lattice 
positions or if they have random size, since Jw , Rw and 
Cu> depend on the distance between the adjacent super- 
conductors and on their geometry 

b. For site disorder the properties of the supercon- 
ducting islands, i.e. their resistance and capacitance to 
ground, are random. More simply, in a dilute array cer- 
tain islands are totally missing. 

Dynamic measurements 10 ' 11 have been done on di- 
lute, or percolative , arrays for p-values somewhat larger 
than the percolation threshold p c = 0.5, where the dis- 
order has a marked influence on the properties of the 
array. At temperatures that are low compared to the 
Berezinskii-Kostcrlitz-Thoulcss (BKT) transition tem- 
perature, Tbkt (which in a dilute array is lower than 
for the regular counterpart, due to the absence of cer- 
tain bonds, see below) the equations of motion can be 
linearized, taking into account only small amplitude ex- 
citations. The problem at hand is then the same as 
determining the vibrational modes of a disordered har- 
monic solid, which is a well studied field. 12 In reference 13 
the impedance Z(uS) of such a dilute JJA has been de- 
termined by evaluating the dynamic voltage correlation 
function by different methods. The coherent potential 
approximation (CPA) has proven to be a useful approach 
for treating dynamic disorder problems. The inverse of 
the impedance is expressed by an effective coupling func- 
tion K(u>,p) depending on frequency and on the dilution 
fraction p : 



Z(w) = fa + 



K(lu, P ). 



(4) 



corresponding to the effective inductance L and resis- 
tance R of the array. 

As a main result for the BKT transition these calcula- 
tions reveal the existence of an effective Joscphson cou- 
pling constant, J e ff- For bonds where one or both of the 
superconducting islands is missing, the coupling energy 
vanishes, i.e. Jw = 0. Thermodynamics is then gov- 
erned by an effective coupling constant </ e // given by the 
relation 13 



j _P_J± 

J eff — ~, J 

1 - Pc 



(5) 



and thus the transition temperature Tbkt {^bTbkt = 
j J e ff or 0.9 J e ff for triangular or tetrangular array re- 
spectively) also decreases. 

Close to and above Tbkt the dynamic response is 
then dominated by the motion of thermally created vor- 
tices. An equation of motion for vortices can be obtained 
starting from the RSJ equations for the superconduct- 
ing phases 3 . Much analytical and numerical efforts have 
been devoted 3 to describing the dynamic behaviour of 



the vortices in regular JJAs. One of the key quantities is 
the frequency and temperature dependent vortex mobil- 
ity /j,(oj,T). For calculating this quantity the Coulomb 
interaction between the vortices has to be taken into ac- 
count, which makes the problem difficult. 

For disordered JJAs this is an even more complex prob- 
lem, since one should now describe the motion of inter- 
acting particles in a random potential landscape. We will 
treat the vortex response mainly above Tbkt where the 
interaction between vortices is screened and will thus be 
neglected. The influence of the random potential on the 
vortex dynamics will be treated in the Multiple trapping 
model, that has been developed for handling the motion 
of electrons in disordered semiconductors. This model 
will be presented in section II. It allows to calculate the 
average vortex mobility /j,(u),T) in the dilute array as a 
function of frequency u> and temperature T. In section 
III the measured electrodynamic response, expressed by 
the complex impedance of the array, will be related to /i. 
In section IV our results for (i will be presented and anal- 
ysed. The main features of our theoretical results for the 
complex impedance of the array are in good agreement 
with experimental data 11 . We shall also show theoreti- 
cally calculated flux noise spectra for disordered arrays. 
Finally in chapter V we shall give some conclusions. 



II. THE MULTIPLE TRAPPING MODEL 
A. The multiple trapping equations 

In this chapter we discuss our model for the vortex 
dynamics in dilute 2d JJAs (missing sites or bonds) . De- 
viations from a regular lattice structure will have an in- 
fluence on the vortices of an array through the Peierls 
force or pinning potential. The equilibrium arrangement 
of the (thermal or field induced) vortices and antivor- 
tices will correspond to a minimum of the free energy in 
the given random pinning potential landscape. By ther- 
mal excitation, in particular above Tbkt, vortices will 
move around in the dilute array thus giving rise to inter- 
esting power law behaviour in the dynamic impedance 
Z(uj). 14 - 16 

Our main goal is therefore to calculate the frequency 
dependent resistance (R) and inductance (L) of a dilute 
array. In calculating R and L we have to first find the 
mobility (/i) of the vortex (V) and antivortex (A) in the 
disordered array. For this we use the so called multiple 
trapping model (MTM), developed for electronic trans- 
port in amorphous semiconductors. 17 In this model the 
regions where sites of the array are missing are regrouped 
into holes of different shapes and sizes. The motion of a 
given particle (we will not distinguish between V and A, 
since in the absence of interaction they undergo the same 
influence by the random potential) is then described in 
probabilistic terms. At a given position r in the array 
one can either be in one of the holes of the array or in 
a regular region, occupied by superconducting islands. 



Thus the state of a particle sitting at r is determined by 
the following probabilities: 

p(r, t) = probability that a given particle is free at 
position r and time t, i.e. it is sitting between holes, in 
a regular region. 

p n {t) = probability that the particle is trapped in one 
of the holes which is indexed by n. 

In the multiple trapping model the time evolution of 
these probabilities is governed by two rate equations: 



dpjr, t) 

dt 



E 



dt 



Vj 



(6) 



with 



dt 



-lr,nPn{t) + lt,nP{*, t) (7) 



Here jt,n is the probability per unit time (transition 
rate) for a particle to get trapped inside the hole n and 
7 ri „ is the release rate out of the hole n. 

If the particle remains in a regular area it contributes 
to the current j = vp(r, t) with its velocity v = /ioE 
determined by the mobility [Iq for a regular array of the 
same structure, and the effective field E which will be 
assumed here as unidirectional (in X direction). Thus 
the first of the two rate equations can be written as 

dp(x,t) j-,dp(x,t) m-^ ^ 



Of 



dx 



Now we use Laplace transformation 

f(z)=^dte^f{t) 

f(0) being the initial condition. 
For (8) this yields 



(8) 



(9) 



zp(z) - p(0) + [i E 



dp(z) 
Ox 



^2 lr,nPn{z) ~ ^ "ft,nf>(z) 
n n 

(10) 



and for (7) 

Zp n (z) = lt,nP{z) - 7r,nPn(z) +Pn(0) (11) 

Imposing the initial condition p(0) = 1 and p n (0) = 
we get from the previous equation 



Pn{z) 



lt,n 



Z + 7r,n 



p{z) 



Now we define 



n n {z) 



It,' 



z + 7, v 



and equation (10) becomes 

zp(z)(l + ir(z)) = 1 - noE 



dp{z) 
dx 



(12) 



(13) 



(14) 



The quantity 



7F ( Z ) =^2 7T n(z) 



(15) 



will turn out to play a central part for the frequency and 
temperature dependence of the electrodynamic response 
of the array. 

In order to get concrete expressions for our transition 
rates we consider thermal equilibrium in zero external 
field where the probabilities, p{v, t) = po an d Pn{t) = Pm 
are independent of space and time. The former is simply 
proportional to the total regular area p. (La) 2 , a being 
the lattice constant, L 2 the total number of sites : 



Po = A.p.(Laf 



(16) 



The probabilities p n are given by a Boltzmann factor, 
involving the binding energy E n of a particle trapped in 
hole number n, and the surface S n of the corresponding 

hole 

Pn = AS n e (iE - (17) 

The thermally activated release rate is given by 

lr,n = re- pE - (18) 

with r being the attempt frequency and (3 — -jr^p- 

In the case of detailed balance all terms in equation 
(7) vanish, therefore 



P0lt,n = Pnlr,n 



(19) 



and 7*,„ = ^ r ,n = r^p- 

For the circular holes, considered here, the area is given 
by S n — S(N) — Na 2 , N being the number of missing 
sites. 

Combining these relations we find the following expres- 
sion for the trapping rate jt.n — lt{N) 



N N 

lt{N) = r— - = r— — 
pL A pN to 



(20) 



B. The vortex mobility 

The mean velocity v of a V is related to its mean po- 
sition x by 



« = 5= id 



d P (r, t) 
dt J " ' ~~dT~ 

Through Laplace transformation (9) we get 
v = zx — / d 2 rxzp(r, z) 



(21) 



(22) 



taking, for simplicity, x(0) = 0, i.e. the vortex started 
(t = 0) its motion from the origin of the lattice. Now us- 
ing equation (14) the previous equation, through partial 
integration, becomes 



d rx 



1 _ /i0 £M£^i 
1 + ttO) 



dx 



l + 7T(z) 

VqE 
1 + ttO) 



2 dp(r, z) 

a rx- 



dx 



d 2 r 



l+7T(z) 



(23) 



The last term of the first part in the previous equation 
represents the initial condition of the velocity and con- 
tributes zero by our assumption. 

The effective mobility /i of a V is given by v = [xE. We 
get from the previous equation 18 



\x(z) 



Ho 



l+7r(z) 



(24) 



a constant r may be understood as arising from the fact 
that our system is considered to be overdamped; hence, 
kinetic energy - which is indeed proportional to tempera- 
ture - is an ill-defined quantity and one may consequently 
be led to assume, that r does not depend on temperature 
and simply represents a characteristic frequency of the 
vortex system. We end up with 



n(u>) 



N„ 



dN- 



D(N)N 



L W) 



o -0E(N) 



(25b) 



The value of N max has to guarantee the correct total 
number of missing sites: 



AS, 



N, 



t (\-p)= J2 D ( N ) N 



(26) 



Therefore N max depends on the value of p, the sample 
size Ntot aud of course on the choice of a hole distribution 
function D(N). 



and tt(z) is given by equation (15). The purpose of our 
calculations is to exhibit, through the mobility /j,, the in- 
fluence of disorder on the dynamic response of the vortex 
system. In order not to mix these effects with the intrin- 
sic frequency dependence of /i in regular arrays and for 
computational simplicity we take /io, the vortex mobility 
in the regular array, to be constant, giving the scale unit 
for h (see expression (38)). 

Taking D(N) to be the number of holes having TV miss- 
ing sites and considering (20) for j t (N), the key quantity 
7r(z = i£l), for real frequencies $7, can be expressed by 

tt(z = m) = w(fi) -» J2n d(n)tt(z = in, n) 



C. Binding energy for circular holes 

The energy of a phase configuration is expressed by 
(2). In the case of percolative arrays 

T _ \ J with probability p , _^ 

"' ~ | with probability 1 — p *• ' 

where J is the Josephson coupling constant between two 
islands. 

Now when phase field 9{r) varies slowly we can write 

1 — cos(8i — 6ii) w v ' 2 1 '' , so we get, converting sum- 
mation into integration for slowly varying phase fields, 



E N Jt(N) 






(25a) 



= K 



A„ 



, v D(N)N 



is a normalized 



In the above expression, D(N) = ^ ' 
number density of holes. We introduce the dimensionless 
frequency lo — -j-, with uj a = kBTB j T ^ a being a char- 
acteristic frequency which includes the lattice constant a 
of the array and the bare vortex mobility /io of the or- 
dered array. Introducing the vortex diffusion constant by 
Dq = ksTBKTlJ'O, for temperature Tbkt, the frequency 
uj a is given by to a = —$■■ Its inverse is thus the average 
time for a vortex to diffuse across one lattice constant 
a. A very sensitive step in our entire procedure is repre- 
sented by the choice of the attempt frequency r. We will 
write the latter as r = u) a f(t) where t = t^- — represents 
the scaled temperature. Here we consider either f(t) — t, 
if the attempt frequency is considered to increase with 
temperature - i.e. r = oj a t -, or f(t) = 1, in the case 
where r is assumed to be constant. At first glance, the 
former case probably appears to be more sound and in- 
tuitive from a physical point of view. However, choosing 



E^Jd z r 



,2 (W) 2 



(28) 



For a vortex at the origin, the phase field is given by 
6{r) = arc tan ^ where x and y are cartesian coordinates 
in the array. We get |V#| 2 = \ and the total energy of 
a vortex in the array is 



J 



1 



E v = - d z r^r 



(29) 



where L is the size of array and n is a lower cut-off cor- 
responding to lattice constant. Now the binding energy 
of the V or A inside a hole (approximated to be circular) 
of radius r^ with -nr\ = (N + l)irrl (the area of the hole 
is equal to N + 1 times the area of a unit cell) is 



E(N) 



J 

J 

' 2 

■■ 7rJln 



1 



d 2 r- 



drr— 

r-2 



J 

2 

2w 



d l r- 



dd-- 

2 



1 



drr— 



2tt 



TrJlnVN + l 



n 



dO 



(30) 



because 7r(r| — r 2 ) = Nnr 2 or — = yW + 1. 
Using (5) 



k B T BKT = -J eff - T -^-J 



(31) 



the exponent of the Boltzmann factor in (25) can be 
rewritten as 



f3E{N) 



1- Pc 2 



P-Pct 

Thus we obtain from equation (25) 
1 



In y/N + 1 



(32) 



tt(u) = 



^dN- ND{N) 



NtotpJi "if + (1 + JV)-°W 

with a temperature dependent exponent 

a(t) = -- 



(33) 



tp-Pc 



III. ELECTRODYNAMIC RESPONSE OF THE 
ARRAY 

A. Resistance and Inductance 

Generally, the response of the array to an external 
electromagnetic excitation can be characterized taking 
the contributions of the superfluid, normal electrons and 
vortices into account. This corresponds to a two- fluid 
model, where the medium is described by an inductive 
superfluid channel in parallel with a dissipativc channel. 



(a) 



R s 



(b) 



Rr 



FIG. 1: Circuit diagram for a 2d superconductor, (a) 
in terms of an inductive and resistive component, (b) 
in the presence of vortices. 



The measured quantity in the array is the sheet con- 
ductance G. The effect due to the vortices can be in- 
corporated in this conductance of the array In the pres- 
ence of a current, the vortices experience a Lorentz force 
which will set them in motion perpendicular to the cur- 
rent flow. Associated with the vortex motion, there is an 
electric field which adds to the electric field of the acceler- 
ated superfluid background. This phenomenon therefore 
leads to an increase in the sheet impedance by an amount 
Zy , which comes in series with the impedance of the su- 
perfluid background, as shown in figures 1(a) and 1(b) 
schematically. 

G is the sheet conductance which is the inverse of the 
sheet impedance Z 



G 



1 


1 


Rs 


iQL s 


1 

Rq 


~T~ {-* SUp 



z- 



(35) 



(34) where G sup is the superconducting part of the conduc- 
tance, the resistance Rq is due to the dissipative processes 
(Ohmic contribution of the conductance) resulting from 
the currents flowing in the junction in the absence of vor- 
tices and antivortices, whereas L s is the sheet inductance 
and R s is the sheet resistance when there are vortices and 
antivortices present in the array. 

We shall now relate the resistive (R s ) and inductive 
(L s ) parts of the sheet conductance G with the resistive 
(Ry) and inductive (Ly) parts of the vortex impedance 
Zy through the following relations : Zy = Z xnn — i$lLQ = 

so 



'jy through the following relations : zjy — ^ sup 



G s up — in,L and G 
Zy = 



__ (-i l_ _1 l_ 

^ Ro R s Rq 



R. 



-r 



R, 



— 



(iOL.. 



— - - ifiLo = Rv 



flL s 

iflLy which 



gives 



R\ 



Ly = 



R s 1 - Rq 

{R7 1 - Rq 1 ) 2 + m s 



1 {VLs)- 1 

CI (R7 1 - Ro 1 ) 2 + (ilL,)- 



-Ln 



(36) 



(37) 



Here Lq, the inductive part, arises from the currents flow- 
ing in the junction in absence of vortices and antivortices. 
The frequency dependent dielectric function of the vor- 
tex system ey(Q) = 1 — iZy/(flL ) is related to the vor- 
tex mobility /i(fi) by the following expression 



ittLoey(tt) = «i!Lo(l + 2irq nii(£l) 



1, 

io,' 



— iftLo + Lo2irq nno( l/ ' '(^) + iv" {Vt)) 
= mL + Zy = iflL + Ry + iflL v {38) 

where go (<7o = 27r J) is the charge of a vortex, n is the 
vortex density and v(fl) = ^-L = v'(Q) + ii/"(Q) is our 
dimcnsionless mobility. 

Solving for the real and imaginary parts from the pre- 
vious equation through the use of equation (5) for the 



triangular array we get 

Ry = 2-Kq 2 jip:§v' L§ = 87Tco a Lo 



P-Pc 



nv'(u>) (39) 



v" 1 — p v" (lS) 

Ly = 2irq 2 nuQ — L a = 87rL -n (40) 

ll p — p c UJ 

where we have used n — na 2 which is the density of V 
per unit cell. 

As Z v + iQL = iflL + 2Trq 2 nL fi (v' '(ft) + iv"(Q)) 
we get 



The temperature dependent expression for Lq, the 
sheet inductance of a single junction, is given by (for 
details see ref. 11 ) 



_i_ = vl| t -<«i-^, 



cVT 



(47) 



Here £ is a critical exponent in two dimensions, b~^ is 
some constant deduced experimentally by setting the en- 
ergy barrier for vortex motion to its theoretical value 10 , 
I c is the critical current of a single junction, c is some 
constant in unit of K -1 ' 2 and T° is the transition tem- 
perature for the superconducting islands. 



1 



1 



1 



1 



R s Rq iVlL s Zy + iQLo 
1 1 



L i(Q + 2irq 2 niJ,ov"(ty) + 27rq 2 nnov' (fl) 
1 a' v -i{Q + a v ) 



L (n + a' v ) 2 +a' 2 



(41) 



where we have used 



a' v = 2irq 2 .nfj,o l/ ' (Q) = (2ir) 2 Jn/iov' (£1) 
o -1-Pc ,, 



P-Pc 



-1/(0) 



2nq 2 n^ v" \Q) = (2tt) 2 Jn^v"{VL) 



— 8 r Kuj a n 



P-Pc 



(42) 



(43) 



where a' v and u' v are the real and imaginary parts of 
the vortex conductance ay which is related to vortex 
dielectric function through ey(Q) = 1 + av ^ - . 

Separating the real and imaginary parts from both 
sides of equation (41) we get 18 



1 

R~ s 

8tt 



1 



R L ((n + a^) 2 + a' 2 ) 

P—Pc V ' 



L u; a ( W + 87T^nu"(uj)) 2 + (87r^nu'(w)) 2 



(44) 



n(n + cr v ) 



L s L (ft + a v f 



V P~Pc V !l 



L Q (« + 87r±5f;m/'H)2 + (Sn^w'(w)) 2 
The scale quantity Rq is expressed by 
Ro = 4.5Rj 



(45) 



(46) 



where Rj is the junction resistance and the prefactor 
4.5 is an estimate deduced by setting the energy barrier 
for vortex motion to its theoretical value (for details see 
ref. 11 ). 



B. Flux noise 

Flux noise measurements give interesting information 
about time correlations in the vortex dynamics. The 
Fourier transform S$(u)) of the dynamic correlation func- 
tion of the magnetic flux threading through a closed loop 
above the array is given by 3 

S*{u) = So J dk y iH : Afc)2 Re \<P PP (k, -iu)] (48) 

where J\(x) is the first order Bessel function and A rep- 
resents the magnetic penetration depth of the JJA and 
4>(k, z) is the Fourier-Laplace transform of the dynamic 
correlator of the vortex charge density py (r, t) 

<t>(k,z)= I ' d 2 r f dte- zt e^- r (py(r,t)pv(O,0)) (49) 

It can be evaluated, for example, by Mori's procedure for 
calculating dynamic correlation functions, which yields 
the following form 



4>{k,z) = 



S{k) 



k B Tk 2 fi(z) 



(50) 



S(fe) 



involving the static charge structure factor S(k) and the 
dynamic vortex mobility p.(z). Neglecting again the effect 
of vortex interaction we use the mobility resulting from 
the multiple trapping model for the flux noise calculation. 
The structure factor should, in principle, be calculated 
by taking into account the effect of the random potential 
landscape due to the holes on the vortex configuration. 
We instead take the simple form 19 



S(k) = 



k Q 



with 



k 2 



2irq 2 n 
k B T 



(51) 



(52) 



which has the correct behaviour for k — > (charge neu- 
trality of the V-A-system) and for k — > oo . 



IV. RESULTS 



In this chapter we present our theoretical results for 
the frequency dependent mobility of vortices or antivor- 
tices in a dilute superconducting array using the multi- 
ple trapping model, the conductance of the array and the 
vortex resistance. Affolter n has measured the properties 
of such a disordered array near the critical percolation 
limit p c below which the superconductivity of the array 
destroys (a publication of his results is in preparation) . 

In order to compare our theoretical results with these 
experimental data we have to determine the values of the 
parameters characterizing the array. 

The measurements n have been performed on a trian- 
gular JJA of lead islands at p — 0.515 with the following 
characteristics: Rj w 5.2777,0, ( ps 1, a s=a 15 x 10 -6 m, 
7 C (0) w 14mA, c w SK 1 ' 2 , T c ° = IK, T BKT = 3.7K, 
b rj 0.22. From this we find: D rj 2 x 10~ 5 ^ and 
uj a rj 10 6 Hz. 

An important goal consists in elucidating the way in 
which the disorder of the dilute array influences vortex 
motion. The key quantity for this is tt(u>) given by ex- 
pression (33). Given the limits of the integral over N 
three frequency regimes can be distinguished, namely 



2-s 



(58) 



(a) w < wi = (1 + N max )- a ^ 
(b) lo>lo 2 = 2- a W 

(c) UJ\ < UJ < Ll> 2 



(53) 



In regions a and b the frequency dependence of the real 
and imaginary parts, 7r'(u>) and tt"{lo), does not depend 
on the details of the hole distribution D(N). It is given 

by 



Region a : 7r'(w) = ttq, 7r"(w) OC lo 
Region b : 7r'(w) oc — j-, 7t"(u>) cx — 



(54) 



In the intermediate region, however, the form of D(N) 
will be reflected in temperature dependent frequency ex- 
ponents for tt'(u) and tt"(lo). A simple estimate can be 
made by assuming a power law distribution of hole sizes: 



D(N) = D N~ 



(55) 



which turns out to correspond quantitatively with the 
arrays, studied in Ref. 11 , where s rj 1.8. For frequencies 
lying well inside the interval [a>i,u>2] the power law of tt 
can be simply calculated: 



rN max ND(N) 

<") a / dN ■ , ,, ,\n- n m ~ 7 °M 



+ (i + iv)- a (*) 



-«(*) 



/o =/ " y TT^ 



(56) 



(57) 
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a(t) 


Wl 


CJ 2 


u(t)(s= 1.8) 


u(t)(s = 0.3) 


0.9 


37 


10 -53 


7.10" 12 


1 


1.05 


2 


17 


9.10" 25 


8.10" 6 


1.01 


1.1 


6 


5.5 


2.10" 8 


0.02 


1.04 


1.3 


10 


3.3 


2.10" 5 


0.1 


1.06 


1.5 



TABLE I: Some parameters for p ■ 

(l + l)" a . 



0.515 and N max = 25; 
= (1 + 25)~ a and uj 2 = 



We summarize in Table I the values of u>\, uj 2 and u 
using the parameters corresponding to the experimental 
data of Ref.ll, namely p = 0.515, 0.9 < t < 10 and 
N m ax — 25. For the given parameter values it turns 
out that the exponent u rj 1 is almost T-independent. 
The resulting frequency dependence the real part of it, 
tt' ~ 1/w, is unusual, since one expects tt' to be an even 
function of u>, as it is indeed the case in regions a and 
b, see relations (54). The power law behaviour (tt ~ 
lj~ u ) is the result of the integration over a contiuous 
spectrum of relaxation rates in (56). For comparison we 
give another form of D(N) which gives more weight to 
large holes (s = 0.3). The corresponding exponent has a 
more pronounced T-dependence. 

The vortex mobility resulting from 7r through (24) will 
determine the electrodynamic response in (44) and (45), 
as well as flux noise (48). Its behaviour, in units of the 
free vortex mobility fio, is shown in figure 2 for different 
temperatures and for different defect concentrations 1 —p. 
The three frequency regimes are again visible, although 
the curves have more structure than the ones for tt since 
tt' {lo) and tt" {m) are combined with each other in the 
expression (24) for v'(io) and v" (uS). 

At very high frequencies, namely u> > 1 the real part 
of the mobility is constant with the bare value [Xq. In 
this regime we are investigating the short-time response 
of the array: the vortices lying outside the holes do not 
have a significant probability of being trapped and their 
mobility is just the bare one. At the other extremity 
of the frequency spectrum, the vortex mobility is also 
constant but the bare value fio is renormalizcd by 1 + ttq 
leading to a reduction by several orders of magnitudes, 
strongly depending on p. Over very long time spans the 
vortices get trapped and released many times and, as 
they do not contribute to mobility as long as they stay 
in a hole, the overall mobility dramatically diminishes 
upon increasing the number and size of the holes. 

The imaginary part of the vortex mobility displays a 
maximum near u> = 1. For lower frequencies Im[n(uj)\ oc 
lo, although at p = 0.6 and below Tbkt a small kink 
appears around lo ~ 10~ 7 (especially for s — 0.3). 

The intermediate-frequency behaviour of the real part 
of the vortex mobility can be understood by looking care- 
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FIG. 2: (ivsu for (a) p = 0.515 (N max = 25) and for 
(b) p = 0.6 (N m ax = 22). s = 1.8. In order to see the 
frequency dependence easily the slopes for u), respec- 
tively lu' 2 , are given. 



fully at the interplay of the real and imaginary parts of 

tt(w) going into /*'(«) = (1+7r< (g'ff „ ( ^ )a ■ 

Two frequency domains show up: In a range below 
u> = 1, extending as far down as 10 -7 for t = 0.9, 
7r"(w) oc lu -1 for all interesting frequencies, while tt'(lu) 



becomes significantly lower than unity. Thus the real 
part of the mobility increases quadratically with fre- 
quency, which is in some sense a surprisingly normal 
behaviour. Conversely, at smaller frequencies, n'(uo) be- 
comes constant, while tt"(u>) oc lu (i.e. it reaches a max- 
imum at some frequency, below the quadratic-behaviour 
window), giving rise to a constant vortex mobility with 
a value intermediate between fig and the low frequency 
limit. The two frequency regimes get shifted towards 
lower w-values when temperature is reduced, as it can be 
seen in figure 2. The curves in figure 2 have been calcu- 
lated for the choice f(t) — t in equation (25b). We have 
verified that the result for /(£) = 1 is almost identical. 

In addition to the intermediate-frequency regime, 
where the real part of the mobility becomes constant, we 
see on Fig. 2, that at t = 2 yet another anomalous regime 
shows up at very low frequencies, where iuf(w) oc u, ap- 
proximately. At high temperature (t — 6) and at fre- 
quencies below the tiny /i oc lu 2 window centered around 
lu ss 1, the mobility plateau gets replaced by a regime 
where fx'(u>) oc w 08 , down to lu *=s 10 -6 . The non-integer 
exponent arises from the fact that, in accordance with 
the Table I, 7r"(w) oc a; -1 ' 04 in this frequency window, 
while tt'(lu) oc w" 1 and tt'(lu) >> 1 still. Obviously, this 
non-integer power-law behaviour of the real part of the 
mobility also arises at lower temperatures, but it becomes 
less pronounced as the exponent u(T) is extremely close 
to unity at t = 2. In order to confirm this behaviour, we 
have computed the vortex mobility at s = 0.3, a choice 
for which the exponent u(T) has a much more marked 
temperature-dependence than at s = 1.8. see Table I. 
In the former case, we find that, at t = 0.9, no relevant 
changes occur, with respect to the results obtained at 
s = 1.8. Nevertheless, at higher temperature t = 2, the 
intermediate-frequency plateau discussed above gets sub- 
stantially reduced and covers at most only one frequency 
decade. At even higher temperature t — 6, the plateau 
is replaced by an anomalous regime extending roughly 
from lu w 10~ 7 tows 10~ 2 , in which the real part of the 
mobility behaves as a power-law lu oc lu~°' 7 again with 
a non integer exponent. For both choices of the param- 
eter s, the real part of the mobility eventually turns to 
a constant at extremely low frequencies. Therefore, in 
conclusion, for temperatures above Tbkt, the mobility 
displays four distinct frequency regimes, which collapse 
into three regimes below the transition temperature. 

We now turn to the electrodynamic response for which 
experimental data are available in ref. 11 . For the vor- 
tex density n showing up in (44) and (45) we have used 
the values obtained in the references 20 ' 21 through Monte 
Carlo simulations of the regular arrays. In order to have 
a consistent treatment one should, of course, know the 
vortex density of a dilute array. The presence of holes 
indeed makes the formation of phase singularities more 
easy, since, in particular, vortices centered in a hole cost 
less energy than in a regular array. However, for tem- 
peratures above Tbkt, where even in a regular array the 
number of vortices grows rather rapidly, the difference 
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FIG. 3: i? s vs u for p = 0.515. (a) s = 1.8 and (b) s = 0.3. 



should not be too important. Our sheet resistance curves 
for p = 0.515 are shown in figures 3(a) and 3(b) for two 
different exponents s in D(N), see equation (55). For the 
sake of comparison figure 4 shows the same curves for a 
slightly larger value of p. The following observations can 
be made: 

- The same three frequency regimes determining the 
quantity 7r(w) and the vortex mobility can be identi- 
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FIG. 4: R s vs lo plot for p = 0.53 (N max = 25) and s = 1.8. 



tied in R s (w) : at very low and at very high frequen- 
cies R s is constant, whereas in the intermediate regime 



it increases as a power of w: R s 



j x( - T \ This tem- 



perature dependent exponent a; is a fine detail of the 
MT model which is in good agreement with the exper- 
imental data. For higher temperature x = 1, which is 
again an interesting signature of the disorder of the ar- 
ray. For lower T the sheet resistance undergoes an up- 
turn to w 2 -behaviour, which is clearly visible in Fig. 3. 
The latter behaviour is intimately linked to the real part 
of the vortex mobility (discussed above) going through 
an intermediate plateau in some frequency window. In- 
deed, within the latter, we notice that Py(uj) <C fjfy(u>), 
trivially implying cr' v (u}) <C ffy(w). Upon considering 
Eqns. (42) and (43) for p close to p c , we furthermore de- 
duce that cr'y(io) is larger than frequency in the "plateau" 
w-window. Consequently, inserting this into expression 



(44), we obtain Rj 1 



Looy(u)' +B.oa' v (uj) 



. Hence, bearing 



R.oLaa v (cj) 2 

in mind that for all the frequencies that we are consid- 
ering here, p! v {uj) oc u>, there exists a frequency window 
in which fj/y(w) = const, implies R s oc cj 2 , especially for 
temperatures close to and above Tbkt, where Lq(T) is 
small. In turn, at higher frequencies, we fall in the u>- 
window where n'y(u}) «w 2 , in which the sheet resistance 
turns to a constant. In conclusion, the anomalous vortex 
mobility plateau which arises from the balance between 
the real and imaginary parts of 7r(w) is at the root of the 
w 2 -upturn of the sheet resistance in our model. 

- In experiment R s approaches a square root frequency 
dependence for u> larger than uj c (see ref. 11 ), rather than 



10 



io v 



x io 6 



J 7 
10' 



/ 




- / y 

/ / 


/■ 




(a) 




t=1.2 
t=2 

— - t=4 
t=6 
t=8 

- - t=10 



10" 



10"-' 10" 

0) 



10'' 



10" 



10' 



.'■ 


,' ' / 

.. / / / 

*' / 
, / / / 

/ / / 


"•" 


- t=1.2 

t=2 

— - t=4 

- - - t=6 
-- t=8 

- - t=10 



10" 



10" 



10" 



FIG. 6: L s vs w plot for p = 0.53. s = 1.8. 
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FIG. 5: L s vs u) plot for p = 0.515. (a) s = 1.8, (b) 
s = 0.3. 



becoming constant. This may point to a new dynamic 
regime that is not fully covered by our model calculations. 
There is however a growing tendency in the Theoretical 
curves to such a further increlase of R s at the highest 
frequencies shown for larger p- values (see figure 4). 

- The temperature variation of the slow frequency level 
of R s , as well as the difference between the low and high 
frequency limits are larger in the MTM results than in 
experiment. This may be a hint that the model, at least 
for the p- value used, attributes too much value to the 
disorder of the array. 

- For higher temperatures the measured data show uni- 
versality: when the curves have reached a linear slope 
they lie on top of each other. In the model results this 
universality is rather well reproduced for a higher value 
of the exponent s of the hole distribution function D(N) 
(figures 3(a) and 4), which corresponds to the effective 
hole hierarchy of the experimental array 11 , giving thus 
less weight to large holes than the lower values of s. 

The inverse of the inductive response, 1/L S , for our 
MTM is presented in figures 5(a) and 5(b) for the same 
parameters as in figure 3 for R s and for a large p in figure 
6. Comparing the two leads to the following observations: 

- The three frequency regimes are again visible: for 
small and for large lo 1/L S varies as uj 2 , whereas in- 
between it is constant. The experimental data 11 show 
this behaviour for low frequency. However, except for 
the lowest temperature, there is no real constant part, 
but relatively smooth cross-over to a square root like be- 
haviour when u> increases. This again points to a high 



frequency dynamics that is not covered by our MT model 
for p — 0.515. However, the experimentally obtained 
y/uj type behaviour is almost achieved in our theoretical 
model for higher values of p and T. 
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FIG. 7: Rv vs l/£ for different Q. (a) for p = 0.515, (b) 
for p = 0.53. s = 1.8. 

On the whole our MTM results are in good agreement 
with the measured data. Certain experimental features 
(in particular the low frequency and low temperature 
behaviour of R s ) are better reproduced for the value 
p = 0.515 of the percolation parameter, which corre- 
sponds to the experimental array. Other details of the 
data are closer to the theoretical curves for p = 0.53. 
This may be due to the oversimplification of the true, 
ramified structure of the regions of missing sites, which 
in the model are approximated by circular hole. On the 
other hand, concerning the probability distribution D(N) 
of these holes the exponent s = 1.8, corresponding to the 



11 



experimental situation, is more satisfactory than a form 
for D which would give more weight to large holes. 

In section IIIA, equations (36) and (37), we have in- 
troduced the notion of a vortex impedance Zy — Ry + 
iflLy. At low frequencies its real part is thermally acti- 
vated as shown in figure 7. This behaviour is confirmed 
by the measurements presented in ref. 11 . The slope of 
the theoretical curves at high enough temperature (e.g., 
around t — 5) yields activation energies of about 4 and 
2.6 in unit of J for p — 0.515 and p = 0.53 respectively, 
in good agreement with the experimental data. For very 
low temperatures the curves become flat. However, in 
between the two limiting regimes there is another com- 
mon slope in the Ry vs 1/i plot at about 1/i = 0.6 for 
p = 0.515 and 1/i = 1.1 for p = 0.53 in our theoretical 
plots. The measurements do not really exhibit our sec- 
ond common slope but the curves in ref. 11 for 317 Hz and 
lower frequencies do indeed show a tendency to a steeper 
slope before they approach the constant value in the low 
temperature limit. In this respect the experimental data 
look more similar to our curves for a larger percolation 
fraction (p — 0.53 in figure 7). 
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FIG. 8: Flux noise S$ vs u for (a) p = 0.515 and (b) 
p = 0.6. s= 1.8. 

Finally we present, in figure 8, the flux noise spectrum 
resulting from our multiple trapping approach, for dif- 
ferent temperatures and p. Flux noise is an interesting 
observable even for regular arrays 3 ' 19 . It is white (fre- 
quency independent) for sufficiently low lu with a value 
that depends on T. For intermediate frequencies the 
curves for different T are identical, having a 1/lu slope. 
This universality has been explained in different ways, 
for example by invoking the dynamic response of bound 



vortex-antivortex pairs which exist below Tbkt, 19 but 
also above, with a finite life-time. It is interesting to 
note that disorder seems to produce a similar behaviour, 
although the intermediate regime has structure super- 
posed on a simple 1/lu dependence, specially for p close 
to p c . In particular, there is another white noise region 
for higher frequency, before the behaviour crosses over to 
1/lu 2 at high frequencies. This structure is again a conse- 
quenc of the frequency dependence of the vortex mobility. 
Moreover, the universality in this regime is not perfect: 
even for p = 0.6 the curves are not on top of each other. 
Unfortunately there are, at least for the time being, no 
experimental data to compare with. 



V. SUMMARY AND CONCLUSIONS 

We have aimed at explaining the electrodynamic re- 
sponse of a Josephson junction array in which a substan- 
tial fraction of the superconducting sites is missing. For 
the motion of thermally excited vortices and antivortices 
in this type of strongly disordered array we have used 
a multiple trapping model. The regions of missing sites 
are grouped into holes which we take to be of circular 
shape. When a vortex or antivortex reaches a hole it has 
some probability of being trapped into the hole (pinning 
effect). It can later on be released due to thermal exci- 
tation. The effect of this disorder is phrased in terms of 
a frequency dependent vortex mobility which determines 
observable quantities, such as the electrical conductance, 
composed by the sheet resistance R s and the inductance 
L s of the array, or the flux noise. We find the following 
results : 

(i) R s exhibits three frequency regimes. At very low 
and high frequency lu a white spectrum is observed with a 
much lower value for low lu than for large u>. For the time 
scale corresponding to the latter regime vortices remain 
either trapped in a hole of free outside any hole. Thus 
the mobility, and the corresponding response is given by 
the one of a regular array. At the opposite end, for very 
low frequencies, vortices get trapped many times during 
one excitation cycle, which strongly reduces their mobil- 
ity. For intermediate frequencies R s oc lu x with x varying 
between x — 1 and x = 2. For higher temperatures the 
R s ex u) almost coincide, signalling some kind of univer- 
sality. These features are in good agreement with the 
experimental data in . 



(ii) For the inductive part we find L s 



CX LU 



2 at low 



lu while for intermediate frequency regime Lj 1 is inde- 
pendent of lu. For large frequencies a tendency of Lj 1 
to grow with frequency with some new exponent is also 
seen for higher values of p. The critical frequency for 
crossover from Lj 1 oc lu 2 to constant Lj 1 decreases with 
the decrease of temperature T and percolation fraction p. 
These results also reproduce well the experimental data 
of reference 11 . 

(iii) For a given frequency the vortex resistance is ther- 
mally activated at sufficiently high temperatures. The 
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activition energy on the order of 4 times the Josephson 
coupling between the existing superconducting sites cor- 
responds to the measured value. For lower temperatures 
a tendency towards another activated form of the Ry 
versus 1/T curves appears, which is also visible as a ten- 
dency in the measured data. 

(iv) We have also evaluated the frequency dependent 
flux noise. l/u> noise is achieved for the intermediate 
frequencies separating the white noise part for small lu 
and l/w 2 noise for very high frequencies, thus the pres- 
ence of disorder leads to similar results as seen in regular 
arrays. 19 . However at very close to p c the flux noise data 
has some unexpected flatcnning before turing to l/u> 2 
part which is a consequence of the frequency dependence 
of the vortex mobility. Unfortunately there are no ex- 
perimental flux noise data for the disordered arrays to 
compare with our theoretical investigations. 
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